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Abstract 

At the critical point in two dimensions, the number of percolation 
clusters of enclosed area greater than A is proportional to A^'^, with 
a proportionality constant C that is universal. We show theoretically 
(based upon Coulomb gas methods), and verify numerically to high 
precision, that C = l/sVSn = 0.022972037.... We also derive, and 
verify to varying precision, the corresponding constant for Ising spin 
clusters, and for Fortuin-Kasteleyn clusters of the Q = 2, 3 and 4-state 
Potts models. 
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1 Introduction. 



It is often useful to characterize critical systems by their geometric properties, 
for example the distribution of cluster sizes which appears to follow a power 
law 

Us ~ Bs-^ (1) 

asymptotically for large s, where gives the number of clusters of s con- 
nected sites, per lattice site. The exponent r is a universal quantity whose 
value is the same for all systems of a given class — for example, ^ for 
all critical percolation systems in two dimensions, no matter what lattice or 
percolation type is considered as long as the rules are sufficiently local. The 
coefficient or amplitude B however is non-universal, varying from lattice to 
lattice. 

Indeed, cannot have a completely universal form because it is written in 
terms of a lattice-level measure, the mass s. Different lattice structures have 
different typical site densities at the lattice level and correspondingly different 
values of B. In order to characterize the size distribution of clusters in a way 
that circumvents the site-level description, the authors of Ref. [|l| considered 
(for the case of two-dimensional percolation) the quantity Nr(£m > i) = N{i) 
which gives the number of clusters whose maximum x- or y- dimension im 
is greater or equal to a given value £, divided by the total system area, 
A = Oi^L"^). They argued that for L >> £ >> a (where a is the lattice 
spacing), this quantity should behave as 

N{£) ~ ^ , (2) 

with the coefficient C being a universal quantity, identical for all 2d perco- 
lating system at the critical point. The universality of C follows heuristically 
from the idea that N represents a macroscopic measure of the large clusters of 
the system, and remains well defined in the limit a — 0, in which the lattice 
disappears. The proportionality to is a consequence of the self-similarity 
of the fractal percolating system, and can also be derived by the following 
argument (in d dimensions): from (|l|) it follows that the number of clusters 
whose mass is greater than s scales as s^""^, and because s ~ where D 
is the fractal dimension of the clusters, the number of clusters whose length 
scale is greater than C scales as £^(1-^), or r'^ by virtue of the hyperscaling 
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relation d/ D = r — 1. This result is valid for any critical system where the hy- 
perscaling relation is valid. Later we shall give other, presumably equivalent, 
theoretical arguments. 

Besides the maximum dimension im, one can consider any other macro- 
scopic measure of the length scale of the cluster, such as the radius of gy- 
ration or the diameter of the covering disk. For each measure, there is a 
corresponding value of C. 

An equivalent way to write (0) is 

N{A) ~ ^ (a^ < A < L^) (3) 

where N[A) is the number of clusters (per unit area) whose area (by some 
measure) is greater or equal to A, and C depends upon the choice of that 
measure. This could be the area of the smallest disk covering the cluster, 
the area enclosed by the cluster, and so on. Eq. (^ is the form of the size 
distribution that will be considered in this paper. 
We note that (D can also be written as P] 

An ~ - (4) 

n 

for 1 << n << where An represents a rank-ordering of the areas, such 

that Ai is the area of the largest cluster, A2 the area of the second-largest, 
etc., for a system whose total area A is defined as unity. Although the rank- 
ordering necessarily starts with clusters whose area is of the order of the area 
of the system, the behavior of @ applies to clusters whose area is much 
smaller than 1 but larger than the lattice element area. Eq. (§) gives the size 
distribution in a proper Zipf 's-law form, in which the weight (here area) is 
inversely proportional to the rank. When written in terms of s, on the other 
hand, the behavior of the size ranking is not a simple inverse power as above 
(compare and |Q) and also is not universal. 

The various measures of the area of the clusters that were considered in 
Ref. included the area of the square Im x Im, the area of a disk that just 
covers the cluster, the area enclosed by the external perimeter (hull) of the 
cluster, and the area enclosed by the Grossman-Aharony (G-A) hull of the 
cluster (in which fjords are excluded) (Percolation hulls are fractal with 
dimension Dh = 7/4 ||^, ^ but enclose a non-fractal, Euclidean area.) For 
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each of these measures, a different value of the constant C apphes, and the 
following values were found: C(square) ~ 0.115, C(disk) ^ 0.104, C(G-A 
hull) ~ 0.037, and C(hull) ~ 0.024. The way these different values of C were 
found was that the first C (square) was measured directly on a fully populated 
lattice (since the measurement of the maximum x- or y-direction is an easy 
task), and then the rest were deduced (in an approximate way) by looking 
at the ratio of the area measures for individually generated clusters. It was 
noticed that C(disk) is close to the fractal co-dimension d — D = 5/48 = 
0.1041666 . . ., but no exact results for any of these C were obtained. 

In the present paper we report on a direct numerical and theoretical study 
of the constant C = C(hull) for the 2d measure of the area enclosed by the 
external perimeter or hull of percolation clusters, Ising spin clusters, and 
Fortuin-Kasteleyn (FK) clusters on the Potts model clusters for Q = 2, 3 
and 4. (Of course, percolation corresponds to the Potts model for Q = 1 and 
the Ising model corresponds to Q = 2.) 

Initially, one of us predicted that for percolation 

C = = 0.022972037 ... (5) 

Independently, the other numerically determined C = 0.022976 ± 0.000005, 
which is completely consistent with this prediction. Additional work de- 
scribed below yields 0.0229723 ±0.0000010 (one standard deviation of error). 
This close agreement confirms that the Coulomb gas methods that are used 
to derive these results are most certainly applicable to percolation and the 
Potts model. We also considered different lattices and types of percolation 
to demonstrate universality. 

For Ising clusters of same-spin sites, we predict the value 

C = — ^ = 0.011486019 . . . (Ising spin clusters) (6) 
16v3vr 

exactly half the value for percolation clusters. For the Potts model with 
Q = 2,3, 4, we also consider the areas enclosed by the FK bond clusters [§], 
and find for the corresponding values of C: 

C = -i- = 0.026525824 ... (FK cluster, Q = 2) (7) 

127r 

a/3 

C = ^ = 0.027566445 ... (FK cluster, Q = 3) (8) 

207r 
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(FK cluster, Q = 4) 



(9) 



47^2 



0.025330296 . . . 



The theoretical justifications of the above predictions are based upon 
considerations reported previously in and expanded upon in the second 
section below. In the third section we describe the numerical work we carried 
out to test these results; we find good numerical confirmation for all the cases. 
Conclusions are given in the fourth section. 

2 Coulomb gas calculation of C. 

In this section we compute the universal amplitude C in the scaling law 
N[A) r^CjA, using Coulomb gas methods [|Ty]. These are not rigorous, but 
are known to give presumably exact results for critical exponents and other 
universal quantities. 

While the development in the Introduction emphasizes clusters and the 
areas enclosed by their external hulls, the focus will shift to both external and 
internal hulls, or loops when both are taken together. A factor of one-half 
will be included in the final results to compensate for this change, so that 
the values for C will be applicable to just external hulls, just internal hulls, 
or the average (but not sum) of the two. 

We consider a finite but large system of linear size L, and total area A = 
0{L^). As will become clear, the precise geometry and boundary conditions 
are not relevant to the calculation of C. L is considered to have dimensions 
of length, so that the total number of sites in the lattice is of order (L/a)^, 
where a is the lattice spacing. 

All the models we consider (percolation, Ising spin clusters, FK clusters) 
are special cases of either the 0(n) model or the Q-state Potts model ||T0(| . 

The 0(n) model is most easily considered on a honeycomb lattice, and it 
is equivalent to the loop gas model defined by the partition function 



where the sum is over all configurations of non-intersecting closed loops on 
the honeycomb lattice. This model has in general two critical points for each 
n in the interval [—2,2]: x = Xci{n) (dilute phase) and x = Xc2{n) (dense 




.total length number of loops 



(10) 



loopconfigs 
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phase). In particular, for n = 1 and x = Xc2 = 1 the loops form the hulls 
of site percolation clusters on the triangular lattice; for n = 1 and x = Xd 
they are the boundaries of critical Ising clusters. For n — > we get a single 
self-avoiding loop, and for n = 2 the loops are the steps on a surface at the 
roughening transition. 

The partition function of the Q-state Potts model is more easily consid- 
ered on the square lattice, and it may be transformed into that of the random 
cluster model, proportional to 

^ ^ ^number of bonds ^number of clusters (H) 

cluster configurations 

where x = p/{l—p). The hulls of the clusters form closed loops on the medial 
lattice, a square lattice whose vertices lie at the midpoint of the links of the 
original lattice. At the critical point x = xdQ) = \jQi Zq is proportional to 
the partition function of a loop gas 

(^g) number of loops (^2) 

loop configurations 

Note that both internal and external cluster hulls are counted as loops. 

In both cases, then, the critical models are equivalent to loop gases with 
fugacity n (resp. \/Q) per loop. The hulls of the FK clusters are in the same 
universality class as the dense phase of the 0{n) model with n = ^Q. 

As discussed in Sec. |I], we are interested in the number N{A) of such 
loops, per unit area of the lattice, whose internal area is greater than a given 
A. Note that we consider A as having dimensions (length)^. For A <^ L^, we 
expect that N{A) has a finite limit as L — oo. In order to obtain universal 
results, we also consider A ^ a^. Our computation of the form of N{A) in 
this regime is in two stages: first we show, from Coulomb gas arguments, that 
the total area contained in all loops behaves logarithmically, oc ^ln(L/a), 
as L/a 00, with a calculable coefficient; then we argue from this that 
in the regime of interest N{A) ~ C/A, with C simply related to the above 
coefficient. 



2.1 Total area inside all loops. 

We shall present two a priori independent arguments, both however based 
on Coulomb gas methods, for evaluating the leading behavior of the the total 
area in side all loops in large but finite system. 
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2.1.1 Wilson loop method. 



The argument of this section follows that of Refs. [0, but, in order to be 
self-contained, we present it again, perhaps with greater clarity. 

The loop gases described above may be mapped exactly onto a height 
model on the dual lattice, as follows. Each loop is assigned a random ori- 
entation, so that a configuration of m unoriented loops corresponds to 2™ 
configurations of oriented loops. There is a 1-1 mapping between the config- 
urations of this oriented loop gas and the heights h, conventionally chosen 
to be integer multiples of vr, as follows: assign h = on the boundary, and 
increase (decrease) h hj n each time a loop is crossed which goes to the left 
(right). The fact that the loops are closed makes this a consistent procedure. 

The weights of n (resp. y/Q) associated with each loop may be taken 
into account in the ensemble of oriented loops in (at least) two ways: the 
most natural would be to assign equal weights |n to each orientation: we 
refer to this as the real ensemble, and denote averages with respect to this 
ensemble with conventional brackets (...). However, these weights have the 
considerable disadvantage of not being local when expressed in terms of the 
height variables. Instead, for calculational purposes, a different weighting 
is usually chosen, in which the phases e^'" are distributed along each loop 
by assigning a phase e^"^/^'^ each time the loop turns leftwards through an 
angle 9. Each anticlockwise (clockwise) loop thus accumulates a total phase 
e'" (resp. e~'"). On summing over orientations, these account for the loop 
weights as long as 

n = y/Q = 2 cos a (13) 

Clearly, these weights cannot now be interpreted as probabilities, and we 
refer to this as the complex ensemble. Averages with respect to this will be 
denoted by [...]. 

Within each ensemble, we may view each oriented loop as carrying a unit 
current in the sense of its orientation. Let J^{x,y) be the corresponding 
current density. For example, for a current directed along a link in the 
positive y-direction, located at x = 0, J^^ = and Jy = 6{x). This current 
density may be used to give a formula for the area of a single closed loop: 

A = —- [ I \x — x'\6{y — y')Jy{x,y)Jy{x' ,y')dxdydx'dy' (14) 



2 

This formula is valid for any non-self-intersecting loop, and is independent of 
its orientation. If however we now consider the same quantity evaluated for a 
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given configuration of a gas of many loops, and we sum over the orientation 
of each loop independently, Jy{r)Jy{r') will average to zero if r and r' are on 
different loops. Thus the expression 

~ll j\^-AKy-y'){Jy{^.y)JyW.y'))d^d\' (is) 

where (. . .) denotes the average over the loop gas ensemble, gives the mean 
total area (Atot) inside all loops. 

What is the current density J^(r) in the height model parametrization of 
the configurations? Let us imagine that the definition of the height function 
h{r) is extended to in such a way that it is constant within each plaquette. 
An obvious candidate is then = {l/TT)e^i^d^h. Within the real ensemble, 
this is clearly correct. It is easy to see that {J^{r)) = on summing over 
orientations of a given loop which passes through r. However, in general 
[J^(r)] 7^ 0. Consider the average over orientations in the complex ensemble 
for a fixed configuration of unoriented loops: 

[J] ocl-e'" + (-!)■ e-'"^0 (16) 

Instead, we have to consider a different operator as representing the current 
in the complex ensemble: oc e^^d^e~'^^"'^^'^ , which now gives 

[j] oc (e-2'° - l)e'" + (e^'" - l)e-'° = (17) 

as required. In Coulomb gas language, j has charge — 2ia/7r. Since there must 
be overall charge neutrality, this is balanced by a charge +2ia/7r distributed 
on the boundary. 

However, the orientation-averaged two-point function {J{r)J{r')) is not 
correctly represented by [j{r)j{r')], since once again this has the wrong 
charge. Instead we should take 

{Ur)Mr)) = X[Ur)Ur')] (18) 

Note that this does vanish when r and r' are on different loops, by virtue of 
(jn\j. The constant A is fixed by requiring that, for a fixed long loop whose 
sides at x and x' are parallel to the y-axis, after summing over orientations, 
{Jy{x, y)Jy{x', y)) = —n6{x)6{x'). The factor n arises from the sum over loop 
orientations in the height model. This gives 

A ((1 - e-^'")e'" - (1 - e2^'^)e-'°) = -n (19) 



8 



so that A = n/(4isina). 

All of this is exact, on the lattice. The weights in the height model 
are local but complicated, involving as they do the phase factors e^'". The 
central assumption of the Coulomb gas approach is that, for the purposes 
of studying the long-distance behavior of correlation functions, they may be 
replaced by the continuum measure exp(— 5), with S = {g/4:'7i) J{dh)^d'\-. 
The parameter g may be determined by a number of methods [l^, [Tl| to 
give g = 1 — a/n, so that n = = —2cos{TTg). The correct branches are 
1 ^ fi' ^ 2 for Xci, and < g < 1 for Xc2- 

Within this free field theory, it is straightforward to compute the cor- 
relation function [Jfi{T)j^{r')] in terms of the Green function G{r — r') = 
[h{r)h{r')] ~ —{1/g) In |r — r'\. First note that 

[^^^)g-2iah(r')A] = f^t^i^^[h{r)h{r'Y] (20) 
p=0 P- 

= f:^-^^^^pG{r-r'Mrr-'] (21) 
p=o P- 

~ (2ia/7rc/)ln|r-r'|[e-2'"'^("')/"] (22) 

= {2ia/7rg)\n\r -r'\ (23) 

The last equality follows because of the way the phase factors enter the sum 
over orientations, so that [e~^'°''(''' ^/'^] = 1. We thus find that 

[Mr)ju{r')] = {2iXa/7T'g)e^,e,„d^d'Jn\r-r'\ (24) 
= (2iAa/7r^,)e,.e.. - ||) (25) 

where we have introduced R = r — r'. After a little algebra we therefore find 

{Ur)Mr')) = k{n) ^^^^-J^'^''' (26) 

This form of the 2-point correlation function of a conserved current is in fact 
dictated by rotational invariance, but it is the coefficient k{n) which is the 
main result: written in terms of g it isQ 

k{n) = — r = cot(7r^) (27) 

7rgsm{7i:g) ng 



^ There is a misprint in the corresponding equation in Ref. 
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Substituting this result into (pTSD, we notice that the result would appear 
to diverge logarithmically at r = r' . However, this is an artifact of the 
continuum approximation: the result ( PBD is valid only for separations |r — 
r'l ^ a. Since the potential divergence is logarithmic, the amplitude of 
the leading term is insensitive to the precise nature of the modification at 
shorter distances, and therefore we may impose a simple cut-off |r— r'| > a on 
the integral. Similarly, the precise form ( p6|) becomes invalid for separations 
0{L) , but the short-distance leading logarithm In a must always appear in the 
form ln(a/L), on dimensional grounds, independent of the precise geometry. 
Thus the mean total area within all loops behaves as 

(Aot) = {k{n)/2)AHL/a) + 0(1) (28) 

where A is the total area of the system. 

As shown in Ref. 0, in any simply connected region TZ the right-hand 
side of (|28|) is proportional to Z]m(l/'^m), where the Am are the eigenval- 
ues of — laplacian in 7?., with Dirichlet boundary conditions. The leading 
term always has the universal logarithmic behavior shown above. Up to a 
non-universal constant which may be absorbed into the cut-off a, the 0(1) 
remainder is universal and depends only on the shape of IZ. For example, 
for a rectangle it is related to modular forms. 

2.1.2 Relation between total area and the mean depth. 

In this section we show how the leading behavior of (Atot) niay also be found 
using methods of conformal field theory on a cylinder, as an extension of the 
results of Ref. |T3[. Let us define, for a given configuration of unoriented 
loops, the 'depth' d{r) of a given site on the dual lattice to be the minimum 
number of loops which must be crossed to connect r to the boundary. That 
is, it is the number of noncontractible loops surrounding r. In the height 
description, it is the supremum of h{r)/TT over all possible orientations of the 
given set of loops. As with h{r), we may extend the domain of definition 
of d{r) to the continuum plane by assuming that it is constant over each 
plaquette of the dual lattice. Then a little thought shows that the total area 
within all loops is simply 

^tot = / d{r)d\ (29) 
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so that we need to evaluate {d{r)). We do this by first evaluating {d{r,r')), 
where d{r, r') is the minimal number of loops which separate distinct points 
r and r' in the infinite plane. By translational invariance this is a function of 
r — r' only. Let us conformally map the plane into a cylinder of perimeter 27r 
via the usual mapping w = \nz. As |r — r'| ^ oo the images of these points 
are far apart along the cylinder, and (i(r, r') is therefore asymptotically the 
same as the number of loops which wrap around the cylinder between these 
points. 

/^From the point of view of the height model with complex weights, these 
loops must in any case be treated separately, since the factors of e^'^"/^^ all 
cancel, so that each orientation is, a priori, counted with weight 1. As is 
well known, this may be compensated for by inserting operators e^'"'^/'^ at 
opposite ends of the cylinder. The free energy per unit length then is — c/12, 
where c = 1 — (6/(7)(a/7r)^. The first term comes from the Casimir effect of 
the fluctuations of the /i-fleld, while the second is the correction due to the 
flux between the charges at either end. This then gives the correct result for 
the central charge c. 

Let us now count the loops which wind around the cylinder with a weight 
n' = 2 cos a', instead of weight n = = 2 cos a. The free energy per unit 
length will be simply modifled by an additional term 

5/ = (l/27r2^)(a''-a^) (30) 

The mean number of loops which wrap around the cylinder, per unit length, 
is then found by taking n'{d/dn') of this expression and setting n' = n. 
Transforming back to the plane, {d{r,r')) is given by this same coefficient, 
multiplying ln(|r — r'\/a). After a little algebra we then find 

{d{r, r')) ~ {k{n)/2) ln(|r - r'\/a), (31) 

where k{n) is given by Eq. (^7^. The logarithmic dependence of this function 



also follows from the work of Ref. |jT^, who derived that dependence from 



general scaling arguments but did not find the coefficient given above. 

This result, of course, applies to the number of loops separating r and r' 
in the infinite plane. However, notice that it has a logarithmically divergent 
dependence on a, which comes from loops which are much smaller in size than 
|r — r'|. This same divergence should arise if we now consider the number 
of loops separating a given point r from the boundary in a large but finite 
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system, for all points whose distance from the boundary is much larger than 
a (but still much less than L). We conclude that 



J {d{r))d\ ~ -{k{n)/2)A\na ~ {k{n) /2)A\n{L/ a), (32) 
where the last statement holds because L is the only dimensionful parameter 



available to compensate a. This gives a second derivation of Eq. (pSl). 

We note in passing that {d{r, r')) = (l/7r)^((/i(r) — /i(r'))^), in the height 
representation. The fact that this behaves logarithmically with |r — r'| is 
consistent with the hypothesis that, in the continuum limit, the heights are 
distributed according to a gaussian ensemble exp ( — (l/27rA;(n)) /(V/i)^rfV) in 
the rea/ ensemble, even though the lattice weights are nonlocal. However, this 
hypothesis is incorrect: as may be shown by extending the above calculation 
on the cylinder to higher moments, the cumulant 

{{h{r) - h{r')f) - 3((/i(r) - h{r')ff ~ const. ln(|r - r'|/a), (33) 

and does not vanish as it would in a gaussian ensemble. However, note that 
this cumulant decays faster than each term on the left hand side, so that 
asymptotically the distribution of d{r) is normal, as was proved by Kesten 



and Zhang |15 . 



2.2 Relation between k and C . 

In this section we first show how the logarithmic behavior of (Aot) provides 
further justification for the assertion that N{A) ~ for <^ A -C L^, 
then show how to relate the coefficients. Recall that N{A) is the number 
of loops with area greater than A, divided by the total area of the system. 
For A ~ L^, this will also depend on L, so let us write it as N{A,L). On 
dimensional grounds it has the form 

N{A, L) = {l/A)F{A/a\ L/a) (34) 

where a is the lattice spacing. For <^ A ^ L^, we expect it to be 
independent of L, but, a priori, it could depend on a. In this regime, let us 
suppose it has the form 

N{A, L) ~ 2C{l/A){A/a^Y (35) 
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where C is a constant and uj is some exponent. This is not of course the 
most general dependence which is possible, but is that which would arise if, 
for some reason, the area scaled non-trivially with a, that is, had a fractal 
structure. Such dependence, with 7^ 0, would for example occur in the 
distribution of masses, rather than of areas, of percolation clusters. The 
form ( p5D should of course connect smoothly onto the behavior for A ~ a^, 
when we expect that — > const., and A ~ L^, where — > 0. 
Now the total area Atot within all loops is related to by 

A,,, = Y^N{A,L) (36) 

A 

Comparing with (pSD, we see that the contribution to the sum from the region 
A L? will exceed 0(ln(L/a) if > 0, and similarly the contribution 
from ^ A <^ l? will violate this bound if < 0. Therefore u; = 0, and 
Ni^-A) ~ 1C I A for <^ A L^. Admittedly, this argument assumes the 
ansatz (|35|), and the reader may be more comfortable with the hyperscaling 
argument put forward in the Introduction. However, independently of the 



validity of (|35|), our argument shows that if N(A) ~ C/A, then the coefficient 
C is related to k{n). For then the leading contribution from the region 
A is 2Cln(^/a^) ~ 4Cln(L/a), so that comparing once again 

with (p8|), 

C = k{n)/8 , (37) 

with k{n) given by (pT]). 

For percolation cluster hulls and FK clusters in the Q-state Potts model, 
we take n = y/Q = —2cosTig in the dense phase < < 1, which yields 
= |, |, |, 1 for Q = 1, 2, 3, and 4, respectively. For critical Ising spin 
clusters, we take n = 1 in the dilute phase where 1 < (7 < 2, so that g = ^■ 
Then by (|27|) we find the values of C given in Sec. |l] and also listed in Table |I| 
(taking the limit in the case Q = 4). The logarithmic corrections that appear 
for the case (5 = 4 are derived in the Appendix. 



3 Numerical results. 

To test these predictions, we carried out numerical studies of percolation on 
square and triangular lattices with both site and bond percolation, and the 
Ising/Potts models on the square lattice. For percolation we considered two 
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cluster type 


C (theoretical) 


C( measured) 


Percolation 
Ising spin 
Ising FK 
Q = 3 Potts FK 
Q = 4 Potts FK 


l/8v^7r = 0.022972037... 
l/16\/37r = 0.011486019. . . 

l/127r = 0.026525824. . . 
y3/207r = 0.027566445 . . . 

l/47r2 = 0.025330296 . . . 


0.0229721(1) 
0.01149(5) 
0.0265 
0.0278 
0.0258 



Table 1: Predicted and measured values of C for various systems. 

ways to generate the clusters: populating the entire lattice, and individual 
hull generation. 

3.1 Bond percolation — Full-lattice population method. 

In the full-lattice population method, we first assign all bonds on the lattice 
as occupied or vacant with probabilities p or 1 — p, respectively, and then 
carry out all possible hull walks around these bonds. These walks go from 
the center to center of each bond along the diagonals, as shown in Fig. 0, and 
turn by an angle -|-7r/2 when the center of an occupied bond is encountered, 
and by — 7r/2 when the center of a vacant bond is encountered. Each walk is 
completed when it returns to its beginning step. 

For bond percolation on the square lattice, we used a square lattice of 
size 512 X 512 with periodic boundary conditions, with p at the threshold 
1/2. We simulated 10^ samples, amounting to a total of 5.2 x 10^^ bonds 



occupied or not. We used the R9689 random number generator of Ref. [|I| 
In the computer program we employed an array of size 2048 x 2048, so that 
we had distinct array locations to represent the bonds and each diagonal leg 
of the hull walks. 

The enclosed area of a hull walk was found "on the fly" by the following 
method: Initially, the area is set equal to zero. When the walk steps to the 
right, the area is increased by one-half the y coordinate at the center of the 
diagonal step (where we had an array point), and decreased by one-half the 
y coordinate when the walk steps left. The zero point of the y coordinate is 
irrelevant, because its value cancels out. The factor of 1/2 comes from the 
fact that each leg of the hull walk changes the x-coordinate by ±1/2; we are 
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Figure 1: Hull paths for bond percolation, with enclosed shaded areas of | 
(top left), 1 (top right), 2 (bottom left) and 2| (bottom right). These are all 
external hulls - the last case also has an internal hull of area | (not shown) . 
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taking the spacing of the bond lattice to be unity. When the walk closes, this 
algorithm gives the area of the enclosed space, with a sign attached: positive 
areas correspond to external hulls that surround clusters, and negative areas 
corresponds to internal hulls (which are of course external to the clusters on 
the dual lattice). 

The smallest area is 1/2; for positive area this corresponds to the hull 
around an isolated site (one with no bonds attached), and negative 1/2 cor- 
responds to the hull inside a square of four occupied bonds, or equivalently 
around an isolated site on the dual lattice. The area of all the hulls are in 
units of 1/2. (Alternately, one could consider the lattice spacing to be ^J2\ 
then the hulls would all have integer areas, and the system area would be 
2L2.) 

Because we use periodic boundary conditions, there is the possibility that 
some hulls could wrap around the torus once or more before closing into a 
loop. The areas for such loops are undefined, unless taken in pairs, but in 
any case we discarded them because we are interested in clusters whose size 
is much smaller than the size of the system. 

We found the statistics for internal and external hulls were identical 
(within numerical error), as one would expect for this self-dual system, and 
took the average of the two. 

For small A we kept track of the quantity A^^ = the number of loops 
(per unit area) whose enclosed area is exactly A, where A = |,1,|,2,.... 
According to (H), this quantity should behave as 

Na = NiA) - NiA + 1/2) ~ ^ (38) 

so that 2A^Na ~ C for large A. The results are given in Table ^ for A < 5. 

To check these results, we derived the exact expressions for A^^ for | < 
y4 < I given in Table ^ These are for an arbitrary bond occupancy of with 
q = 1 — p. For A = | . . . 3 these expressions are identical to the expressions 
for the number of clusters (per site) containing h = 2A — 1 bonds, which are 
well known Il8[. For larger A we had to make modifications to the bond 



cluster expressions to take into account graphs that contain internal open 
spaces with vacant bonds, which result in areas larger than {b + l)/2. We 
subtracted the term 2p^q^^ from Nz and added it to Na to account for the area 

2 

of an open 1x2 rectangle, whose external hull area is 4, not 7/2. Likewise, 
the term 20p'^q^^ (the 1x2 rectangle with an extra bond attached) was 
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Na{p) 

2pq^ 

p^{Aq^ + 18gi°) 
p\q^ + 32g9 + 55giO) 
p5(8gi° + 30gi2 ^ ^ggg^^ ^ 174^14) 
/(12g" + 40g^2 ^ 332^14 + 572^15 ^ 570gi6) 
2p6^ii ^ p7(2gio + 136^13 ^ igggU ^ 33g^i5 ^ 2030^^^ + 2712g^^ 
+1908gis) 

20j97gi3 + p8(22gi2 + I86gi4 + 844gi5 + 8685^*^ + 4064gi^ + 9972^^^ 
+10880gi9 + 6473g20) 



Table 2: Exact results for Na{p) for bond percolation on the square lattice 

1 9 

2 • • • 2' 



at occupancy p = 1 — q, for A — ^ ^ 



subtracted from N4 and added to Ng/2- Finally, the terms A2p^q^'^, llAp^q^^, 
and p^q^^, which correspond to various graphs with area greater than 9/2, 
were subtracted from Ng/2- These various diagrams are shown in Fig. ^ 
This shifting of terms has the effect of making A^^ follow asymptotically the 
exponent —2 of (|38| ) rather than the exponent — r = —2.055 . . . followed by 

Us- 

Taking p = 1/2 and multiplying by 2A'^, we arrive at the estimates for 
C listed in Table ^. The agreement with our numerical results is excellent 
— within the small statistical error. Interestingly, the convergence of these 
estimates is rather quick — already, at A = 5, the result is within 6% of the 
(presumably) exact value. 

To analyze the data for larger A, we considered the quantity N{A, 2A) = 
the number of clusters whose enclosed area is greater or equal to A and less 
than 2 A. According to (j^), this quantity should behave as 

NiA, 2A) = NiA) - Ni2A) ~ ^ (39) 

so that 2AN{A, 2A) ~ C for large A. 

The measured values of 2AN{A, 2A) are given in Table ^. They mono- 
tonically decrease to a value 0.0229860(45) for A = 2048, but then slightly 
increase at A = 4096; for larger A, the increase continues, as seen in Fig. 
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2p6qll 

A = 4 
20p7ql3 

A = 9/2 



42p8ql4 
A = 5 



Il4p8ql5 



A = 5 



p8ql6 



A = 13/2 

Figure 2: Clusters contributing higher-area terms to polynomials in Table ^ 
Solid hues represent occupied bonds, dashed lines are vacant bonds, and the 
dotted lines trace out the external hull. These are the graphs that have to be 
"moved" in the usual cluster polynomials from A = (6+ l)/2 (where b is the 
number of bonds) to higher A due to the existence of enclosed open spaces. 
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0.02315 




Figure 3: Plot of 2AN{A, 2A) vs. ^-o.875 bond percolation on a square 
lattice. Upper data points: lattice population method. Lower points (shifted 
down by 0.00005): single hull generation method. 
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A 


Full 


Single 


Exact 




lattice 


hull 


results 


1/2 


0.0312500(1) 


0.031247(3) 


1/32 = 0.03125 


1 


0.0312500(1) 


0.031252(4) 


1/32 = 0.03125 


3/2 


0.0263674(1) 


0.026363(4) 


27/1024 = 0.026367188 


2 


0.0253907(2) 


0.025398(5) 


13/512 = 0.025390625 


5/2 


0.0257491(2) 


0.025744(6) 


3375/2^7 = 0.025749207 


3 


0.0254749(3) 


0.025477(6) 


3339/2^^ = 0.025474548 


7/2 


0.0249188(3) 


0.024917(7) 


104517/2^2 = 0.024918795 


4 


0.0249898(4) 


0.025004(7) 


6551/218 = 0.024990082 


9/2 


0.0247714(4) 


0.024778(8) 


13298985/2^9 = 0.02477129 


5 


0.0245659(5) 


0.024557(8) 





Table 3: Values of 2A^Na for small A for bond percolation on the square 
lattice: two algorithms and exact results. Errors in last digit are given in 
parentheses. 

^, (upper curve) where the data from A = 128 to 16384 are shown. We 
attribute this increase to interference of clusters with themselves around the 
periodic boundary conditions, and thus ignore these data. Fitting the the 5 
data points from A = 128 to 2048 as a function of A~^ to a straight line, 
we find a good linear fit with 6 = 0.875 as shown in that figure, with the 
equation of the line given by 

2AN{A, 2A) = 0.0229712 + O.OIUSA-^-^'^^ (40) 

implying C = 0.0229712. 

We estimate the error in the above value of C to be ~ 10~^ from the sta- 
tistical error of the data and the uncertainty in the extrapolation to infinity. 
The predicted value (|]) falls within these error bars. 

In terms of the length scale £ ~ Az, this exponent corresponds to a 
correction of the order f"^'^^, which is the scaling of the hull of the cluster. 



Indeed, this finite-size correction can be interpreted as a surface effect 
reflecting the arbitrariness in locating where precisely the hull of the cluster 
should be placed. 
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As a test of our procedure, we also compared our measurement of the 
total number of loops (hulls of both type) per unit twice the number 

of clusters) with the theoretical result, which for bond percolation on the 
square lattice is given by the Temperley-Lieb result pO|, ^ 

J2Na-- 3^3 - 5 = 0.196152422 . . . (41) 

A 

Our measured value was 0.1961572(14), larger than the above prediction 
by only 0.0000048(28). This difference corresponds to an excess number 
of 1.2 loops per lattice (found by multiplying the latter number by 512^), 
which is barely discernible above the statistical error of ±0.7. In fact, this 
correction can also be predicted theoretically. For a system with a rectangular 
boundary of aspect ratio r, the excess number of clusters is a known function 
b{r) ||2l|, To find the excess number of loops, note that the quantity 
Uc + ric' — Til (the number of clusters, plus the number of dual lattice clusters, 
minus the number of loops) equals 1 if there is a cross-configuration on the 
lattice or the dual lattice, and zero otherwise. Thus it follows that the excess 
number of loops is just 26(r) — 27r+(r), where vr+(r) is the cross-configuration 
probability, which has been calculated by Pinson [^. For a square system, 
the excess number of loops is predicted to be 

2[6(1) - 7r+(l)] = 2(0.883576 - 0.309526) = 1.14810 (42) 



using 6(1) from p2| and 7r+(l) from |jl|. This prediction happens to coincide 



almost exactly with the measured value (even though the error bars of the 
latter are quite large). This predicted value can be tested to higher precision 
most easily by going to smaller lattices. 

Besides the problem of clusters interfering with themselves, there is also 
the problem in the population method that the statistics for larger hulls are 
rather poor because of the relatively small number of such hulls that are 
generated. In the next section we consider a method that addresses both of 
these problems. 

3.2 Bond percolation — Single hull generation method. 

It is well known that percolation clusters can be grown individually through 
a process where bonds are made occupied or not only when they are en- 
countered (the "Leath" method). In the same way, percolation hulls can be 
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A 


Pull 


Single 




lattice 


llTlII 


1/2 


0.0625001(1) 


0.0625022(42) 


1 


0.0429689(1) 


0.0429634(32) 


2 


0.0306645(2) 


0.0306677(26) 


4 


0.0270220(2) 


0.0270226(25) 


8 


0.0250527(3) 


0.0250528(25) 


16 


0.0240647(4) 


0.0240635(25) 


32 


0.0235504(6) 


0.0235463(26) 


64 


0.0232785(8) 


0.0232740(27) 


128 


0.0231360(11) 


0.0231412(28) 


256 


0.0230598(16) 


0.0230616(29) 


512 


0.0230218(22) 


0.0230212(31) 


1024 


0.0229968(32) 


0.0229948(32) 


2048 


0.0229860(45) 


0.0229849(33) 


4096 


0.0229882(63) 


0.0229785(35) 



Table 4: Values of 2AN{A, 2A) for bond percolation on the square lattice 
for the two algorithms. 
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generated individually on a blank (undetermined) lattice by a kind of grow- 



ing self- avoiding walk that mimics the walk used to trace out hulls |24]. For 



critical bond percolation the walker moves along the edges of a square 
lattice (the diagonals in Fig. |l[), and turns by +7r/2 or — 7r/2 randomly at 
each vertex, except at vertices previously visited, where it always turns to 
avoid retracing itself. The walk terminates when it returns to the origin and 
cannot proceed further. Note that pseudo-random numbers are generated 
only for the bonds that are visited during the walk, making this method ef- 



ficient. This walk has also been studied as a kinetic Lorentz-gas model 
and the results here apply to that model also. 

In order for the contribution of a given hull to be the same as on the fully 
populated lattice, it is necessary to weight each walk by where t is the 
number of hull steps. This compensates for the fact that a hull of t steps is 
generated with t times the probability in the single cluster method compared 
to the population method, because there are t places a given walk can start 
from. 

This weighting can also be checked as follows: The probability of gener- 
ating a closed hull of at least t steps is generated with a probability 

P{t) ~ cr^/' (43) 

where c is a constant. Defining a Euclidean length scale i ~ t^^^" , it follows 
that the area enclosed by the walk scales as A ~ £^ ~ f^/Dn _ ^8/7 _ ^^j^ug^ 
the probability of growing a walk enclosing at least area A scales as A~^/^, 
so that the probability of growing a walk of exactly area A scales as A~^^^. 
When we weight a hull by the factor 1/t ~ A^'^^^, we thus get the proper 
probability ^-9/8-7/8 = ^-2 gj^g^^ -^^ 

In our simulations, we considered square lattices of size L x L with peri- 
odic b.c. This is the lattice of the hull walks, which is rotated by 7r/4 from 
the square bond lattice, and has a spacing that is V2/2 of the bond lattice 
spacing. Note that the square system boundary here corresponds to a dia- 
mond on the bond lattice. We stopped all walks that did not close by 65536 
steps, and kept track of the areas of all the walks that closed before this 
cutoff, without wrapping around the periodic b.c. With this cutoff, we could 
be assured that all walks that were stopped at the cutoff would ultimately 
enclose an area of at least 65536/8 = 8192, taking into account that there 
are at most 4 hull steps around each wetted site, and each square on the 
hull-walk (rotated) lattice corresponds to an area of 1/2. 
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While the statistics of walks of areas smaller than A = 8192 should thus 
be unbiased by having this cutoff, they can still be biased by the finite-size 
of the lattice. For runs on lattices of size 1024 x 1024 and smaller, we found 
both wraparound clusters and large deviations in the hull statistics for larger 
A. Even for lattices of size 2048 x 2048, where no wraparound occurred with 
this cutoff, we still found significant, obviously finite-size deviations even for 
N{A,2A) for A below A = 8192. We attribute these deviations to hulls 
making contact with themselves around the periodic b.c, without actually 
closing to wrap around. Therefore, to be absolutely certain of no finite-size 
effects, we went to a lattice of size 65536 x 65536 using the virtual lattice 
method of [^. We checked that with the cutoff of 65536 steps, indeed no 
walk got anywhere near the boundary of the system. 

We carried out 1.8 x 10^ walks on this lattice, which, like the simulations 
for the 10'' fully populated lattices, took several weeks of workstation com- 
puter time. A total of 3.2 x 10^^ hull steps simulated here, compared with 
1.0 X 10^^ in the simulations of the populated lattices. The algorithm for the 
single-hull method is somewhat simpler and more efficient than that for the 
lattice population method. 

In the single-hull method, larger hulls are generated with a higher prob- 
ability than in the lattice-population method: the number generated in the 
interval {A, 2A) (before reweighting) is proportional to A~^^^ here, compared 
with A~^. This is advantageous because the large hulls with their small finite- 
size effects are essential for finding C accurately. On the other hand, in the 
single-hull method a large fraction of time is spent on the walks that reach 
the cutoff before they close (and are discarded): Eq. ( ^31) implies that the 
total number of steps for all the hulls that reach the cutoff tmax grows as 
~ ct^l^, while the total number of steps for all the hulls that close before 
tmax is given by 

■'-«(-f)*-^^e. (44) 

Thus, no matter what the value of the cutoff is, a fraction 6/7 = 85.7 % 
of the work (ignoring finite-size effects) is spent generating walks that reach 
the cutoff without closing and are thus discarded. Still, for very large cutoffs 
this overhead is compensated by the increase in useful statistics for large A, 
making this method advantageous. 

Note that, in our simulations of 3.2 x 10^^ hull-walk steps, the fraction 



24 



of those steps belonging to clusters that reached the cutoff tmax = 65536 was 
6.000124/7, with the deviation from 6 in the numerator being about equal 
to the apparent statistical error, ^ 0.0001. This result seems to provide a 
very precise confirmation that the exponent in P{t) is indeed —1/7 (i.e., the 
hull fractal dimension is Dh = 7/4), although to quantify the precision of 
this result one would have to investigate different values of the cutoff tmax to 
determine the finite-size corrections. 

For small A, results for 2A'^Na are given in Table ^ agree with the exact 
values, confirming that the 1/t weighting is correct. Because the single hull 
method gives fewer of these small hulls than the lattice population method, 
these results have larger error bars. Here use used (iV^(total))~^/^, where 
iV^ (total) is the total number of clusters of size A, to estimate the error 
bars. 

Likewise, the results for N{A, 2A) for all A, given in Table ^ are seen 
to agree with the lattice population results. For the largest size ranges, the 
single-hull method is seen to give better error bars (and are not biased by 
finite-size boundary effects). 

A plot of the 2AN{A, 2A) vs. A-^-^'^^ for 128 < A < 4096 is also given in 
Fig. |l| (shifted down by 0.00005), and the data are fit by the linear function 
given y 

2AN{A, 2A) = 0.0229692 + 0.01197^"°-^^^ (45) 

which is consistent with the results of the lattice population method (^). 
The error bars on the intercept is about the same, 10~^. 

Thus, although the single-hull method is in principle advantageous, for 
the system size we considered we obtained C with about the same precision 
as the lattice population method, with about the same amount of work. 
However, the single-hull method allowed us to show that the curvature in 
the behavior of 2AN{A, 2A) for large A as seen in Fig. ^ was indeed due to 
cluster interference around the periodic boundaries. 

Assuming the predicted value of C given by (^, we can also make a plot 
of log(2AN(A, 2 A) — C) vs. \og{A) (not shown); with single-hull data we find 
good linear behavior with a slope —9 = — 0.88±0.01, which is consistent with 
the value 0.875 that we have been using. 
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Figure 4: Medial lattices used for hulls in site percolation: (a) square lattice, 
(b) triangular lattice, and (c) brick-lattice form of triangular lattice. 
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3.3 Site percolation on the square and triangular lat- 
tices. 



We also carried out simulations of site percolation on two different lattices 
to demonstrate the universality of the result for C. 

For site percolation, the logical choice for the hull walk around a cluster 
is to follow a path on the medial lattice whose vertices are at the center 
of the faces of the lattice, as shown in Fig. ^ for the square and triangular 
lattice. This choice allows the single isolated site to have a non-zero area, 
and is symmetric for internal and external hulls for the triangular lattice. 

For the square lattice, we carried out 4 x 10^ samples on a lattice of 
size 256 x 256, using the weighted single- hull method. (In our program we 
employed a computer array of size 512 x 512 to include the sites of the medial 
lattice.) With such a small lattice, finite-size effects appeared for hulls with 
A larger than ^ 1024. We used occupancy probability p = 0.592746, which 
is close to the critical threshold for this system [^]. Here we generated the 



hulls starting from a segment between a single occupied and vacant site, 
which occurs in a populated system with a probability of p{l ~p). The latter 
factor was therefore included in the total weight of each hull, along with the 
1/t weight, where here t is the number of steps along the medial lattice. 

We found that the statistics for internal and external hulls are quite 
different, as one would expect by the asymmetry of this system. For example, 
for A = 1, Ni = p{l — p)'^ ^ 0.0163053 for an external hull, and A^^i = 
(1 — p)p^ ~ 0.00620604 for an internal hull. This large difference persists 
as A increases, and suggests that some other definition of the hull which 
gives more symmetric results between external and internal hulls might be 
advantageous. 

In Fig. ^ we show 2AN{A, 2A) for the two kinds of hulls, along with 
their average. Taking the average is the same as including both types of 
hulls in the area calculation (and dividing by two). Indeed, in the theoretical 
development in Section ^ both internal and external hulls were included 
in the calculation, so it is appropriate to take this average. The finite-size 
corrections to the average measure again followed a behavior with exponent 
close to —0.875, which was used in the plot in Fig. ^ The line in that figure 
is fit by the equation 

2AN{A, 2A) = 0.022976 - 0.0114A"°-^^^ (46) 
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Figure 5: Plot of 2AN{A, 2A) vs. yl^o-^^s f^j. gj^^g percolation clusters on a 
square lattice: external hulls (triangles), average (circles) and internal hulls 
(squares). The equation of the line fit through the average points is given in 
Eq. (il). 

where the unit of area is one square lattice spacing on the square lattice. 
The average measure extrapolates (for large A) to a value ~ 0.022976, in 
obvious agreement with the theoretical prediction, and making a more precise 
determination rather superfluous. 

Note that, the coefficient to the correction term is similar in value to the 
coefficient for bond percolation, even though a different kind of path was used 
to define the hulls in the two cases. The similarity might be a coincidence, 
or it might reflect a fundamental equivalence of perimeter corrections for site 
and bond percolation on this lattice. 

For site percolation on the triangular lattice the medial lattice is a hon- 
eycomb lattice with a hexagon around each vertex of the triangular lattice 
as shown in Fig. |^. To implement this in the computer, we used the square- 
lattice form of the honeycomb lattice, also shown in that figure, where the 
hexagons become rectangular bricks and a single site on the triangular lattice 
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now becomes a pair of sites on the square lattice. Thus, we could use the 
same basic algorithm as we used for site percolation on the square lattice, 
with the only modification being that sites are occupied or made vacant in 
pairs, with a probability 1/2. 

For this case (the triangular lattice) we used the lattice-population method 
on an underlying square lattice of size 1024 x 1024 with periodic b.c. We 
generated 2.4 x 10^ independent samples. As expected, the internal and ex- 
ternal hulls had equal statistics, within error, reflecting the symmetry of this 
system. 

Again, the data closely followed the ^-o-S'^^ ^gi^avior, and we do not plot 
it. Fitting the data in the range 2^ < A < 2^^ (where A is measured in 
square-lattice units, so that the smallest hexagon corresponds to A = 2) we 
found the following behavior: 

2AN{A, 2A) = 0.022977 - 0.0146^-°-^^^ (47) 

again agreeing with the predicted value of C. In this case, that value is 
approached from below for finite systems, with the definition of cluster-hull 
area used here. 



3.4 Ising clusters. 

To study the clusters of the Ising model, we considered a square lattice of 
size 1024 x 1024 with periodic b.c, and simulated the system at the critical 
temperature of exp(— J//cT) = 1 + v^2 (where J is the coupling constant in 
the Potts model formulation, H = —JJ^^ai^aj) using the Wolff variation 
^ of the Swendsen-Wang method 0]. We initialized the lattice with 1000 
updates, and then measured the hull area distribution treating the system 
exactly as if it were one of site percolation, using the same definition of hull 
areas as shown in Fig. 2 and indeed the same algorithm. This was followed by 
10-100 Wolff updates, and the procedure was repeated. 140,000 realizations 
were generated. 

As in the site percolation case, we found rather large differences between 
internal and external hulls, as seen in Fig. |^. Here we also found very large 
deviations for large A, presumably reflecting stronger correlations due to the 
interaction. (Indeed, runs on a smaller 256 x 256 lattice showed even stronger 
large- A deviations.) The average measure for the smaller hulls is consistent 
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0.014 



0.0135 




Figure 6: Plot of 2AN{A, 2A) vs. A~^-^'^^ for Ising clusters: external hulls (tri- 
angles), average (circles) and internal hulls (squares). The line is fit through 
the four rightmost points of the average values. 
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with the A ^ finite-size scahng used in that figure, and a fit of the points 
for A = 2 through 16 yields the straight fine as shown in that figure, with a 
fit of 

2AN{A, 2A) = 0.011487 + 0.004458/1-°-^^^ (48) 

The intercept is nearly identical with the predicted value of C d), in spite 
of the rather small size of clusters that were used. We estimate the error to 
be ±0.00005. 



3.5 FK Clusters on the Potts model for Q = 2, 3 and 
4. 

We also studied the FK clusters on the Potts model at the critical temper- 
ature e"^/*^"^ = 1 + y/Q. These clusters are the bond percolation clusters 
when bonds are drawn between neighboring identical spins with a proba- 
bility 1 - e--^/'^^ = x/Q/(l + ^Q). We defined the hulls exactly as in the 
square-lattice bond percolation case (Fig. 1) and indeed could use the same 
algorithm to trace out and measure the hulls after the bonds have been spec- 
ified. 

To thermalize the system, we used the Swendsen-Wang (SW) procedure 
of identifying all FK clusters on the lattice and then randomly reassigning 
their spins. Indeed, the FK hull measurements and the SW update method 
naturally go hand in had in this calculation, since the identification of the 
FK clusters is needed for the SW method. For Q = 2 and Q = 3 we used 
a lattice of size 512 x 512 and obtained the results shown in Figs. |^ and ^ 
where once again we find large discrepancies between internal and external 
clusters, and take the average of the two. That average is found to fall on 
a nearly straight line when plotted as a function of A~^ taking 6 = 0.875 
for Q = 2 and 6 = 0.7 for Q = 3. The extrapolated exponents are seen to 
approach the expected theoretical values, as shown in Table |l]. We simulated 
82,000 samples (Q = 2) and 1,000,000 samples (Q = 3). 

Note that the discrepancy between internal and external hulls refiects 
an inherent asymmetry for FK clusters of the Potts model in finite periodic 
systems for Q > 1. This asymmetry is also manifested in the behavior of 
the fraction of bonds that are occupied, which in finite systems has a value 
somewhat greater than the infinite-system value of | [31 . 



For (5 = 4, very large differences between internal and external hulls 
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Figure 7: Plot of 2AN{A, 2A) vs. for FK clusters of the Ising model 

(Q=2 Potts): external hulls (triangles), average (circles) and internal hulls 
(squares). 
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Figure 8: Plot of 2AN{A,2A) vs. A'^-^ for FK clusters of the Q=3 
Potts model: external hulls (triangles), average (circles) and internal hulls 
(squares) . 
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persisted even for relatively small values of A on the 512 x 512 lattice, so 
we went to a larger lattice of size 2048 x 2048 (20,000 realizations) which 
improved the behavior somewhat. Even for this lattice, however, large finite- 
size effects were apparent. Similar large finite-size corrections have been 
seen in other Potts model studies at Q = 4 (e.g., |^ ^) and are generally 
expected to be logarithmic in character. In the Appendix we have calculated 
these corrections analytically for this case and find 

iV(A)~^(l-^ + 0((lnAn) (49) 

where 02 is a constant. The above result implies that 2AN{A, 2A) = C + 
0{(\nAy^). In Fig. | we plot our resuhs for 2AN{A,2A) function of 
{lnA)~^. The data fall on a straight line for large A, and the intercept yields 
C = 0.0258, which is comparable to our predicted value of l/47r^ = 0.0253 . . .. 

Note that if we plot the data versus 1/ In A, we find about as good of a fit 
to linear behavior for large A, but then the intercept would 0.0231, quite a 
bit below the predicted value of C. Likewise, if we fit the data to a power-law 
as we did for other values of Q, we find fairly linear behavior with an abcissa 
of A~^-^, but now the intercept is 0.0279. Thus, the data is consistent with 
our prediction for C combined with the predicted l/(ln A)^ finite-size scaling. 



4 Conclusions. 

We derived and numerically confirmed predictions for the behavior of the 
area-size distribution of various Potts model including percolation clusters. 
For the latter, we also considered different lattices and percolation types 
(site and bond) to demonstrate universality. The theoretical ideas presented 
in Section ^ were well verified numerically, especially in the percolation and 
Ising model cases. For Q = 4, our results were consistent with the logarithmic 
finite-size behavior predicted here. 

This work confirms the idea of a universal size distribution expressed by 
Eq. (^). An alternate way to state that result is as follows: Consider that the 
unit of area is now some value A much smaller than the lattice size (which 
is therefore no longer of unit area). Then, (^ implies that the number of 
clusters whose enclosed area is greater than A, per unit area A, is a constant 
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Figure 9: Plot of 2AN{A,2A) vs. (InA)-^ for FK clusters of the Q=4 
Potts model: external hulls (triangles), average (circles) and internal hulls 
(squares). 
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C, for all values of A. The lack of dependence on A is a direct consequence 
of the scale-free nature of this fractal system. 

The arguments put forward in Sec. [^.1.2| also imply that the number of 
cluster hulls which must be crossed to connect a typical point deep inside 
the system to the boundary behaves as 4Cln(L/a), where L is the system 
size, with same value of C for each universality class. So, for example, the 
fact that C for critical Ising spin clusters is half that for percolation clusters 
means, according to Zipf's law, that the nth largest cluster is the Ising case 
has roughly half the area of the nth largest percolation cluster. This is 
consistent with the fact that we have to cross one half as many cluster hulls 
to reach the boundary in the Ising case. It might suggest that we may go 
from the ensemble of percolation hulls to those of Ising clusters simply by 
erasing every other percolation hull, e.g. by ignoring all the internal hulls! 
This however is not the percolation hulls have a different fractal 

dimension from those of Ising clusters. 

The form of is also consistent with the existence of the universal 
amplitude ratio, = [a(l — a)(2 — a)J-'cY^'^$,o, where a is the free-energy 
critical exponent, J^c is the critical part of the free energy per unit area, and 
^0 is the amplitude for the correlation length |3^. For any value of Q in 
the random cluster model, dJ-'/dQ gives the mean total number of clusters 
per unit area J2s^s = Y^a^a- At the critical point, Na ~ —N'{A) ~ C/A"^ 
for yl ^ a^, and near the critical point one expects a scaling law A^^ = 
A~'^^{A/C,^), where is some nontrivial scaling function with $(0) = C, 
which decays exponentially fast as m — > cxo. This gives, on substitution into 
EaNa-^I^NacIA, 

^ n, ~ const. + BC"^ (50) 

s 

where the constant is nonuniversal, as it depends on the details of the cutoff, 
and 

B^r^^iu. (51) 

Eq. ( pOD is of the form expected from hyperscaling [Q, with B = (R^)^ /]a{l— 
a) (2 — a)] directly related to the universal combination (recently given 
exactly for percolation by Seaton However, we see that it is given by a 

certain integral over a nontrivial scaling function, while C is just one limiting 
value of this function. 
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The results presented here represent the first examples where a measure 
of the cluster size distribution is given exactly (in the asymptotic limit), both 
in exponent (here, simply —1) and amplitude (the value C). The agreement 
between the theoretical prediction and the numerical results for percolation 
(to a relative accuracy of better than 10"'*) compares well with other pre- 
cision tests of conformal field theory predictions for percolation amplitudes, 
for example the crossing formula (where the results have been confirmed 
within a relative error of about 10~^ 0). Knowing the exact result for C 
at the critical point allows finite-size effects and behavior away from the crit- 
ical point to be studied, without at the same time having to determine these 
critical parameters. In percolation especially there has been great interest 
in size distributions and their finite size corrections, so this result should be 
useful in that field. 
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Appendix. Logarithmic corrections for Q = 4. 



We summarize the arguments leading to Eq. (^). It has long been known 
that many critical quantities in the 4-state Potts model exhibit confluent 
logarithmic corrections. In the RG framework, this is explained by the exis- 
tence of a marginally irrelevant scaling variable |^|. A general formalism for 
computing the form of these corrections was developed in Ref. was taken 
further in Ref. and recently has been applied to the fractal properties 
of Q = 4 FK clusters by Aharony and Asikainen ^T[|. In general [^], loga- 



rithmic corrections to susceptibilities take the form of multiplicative powers 
of logarithms, and are therefore numerically very significant, but in some 
quantities, for example the finite-size scaling of the free energy at the critical 
point ^2|, they give only additive corrections. We shall argue that this is 



the case here. 
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Following suppose that the fixed-point hamiltonian is deformed by a 
marginal perturbation H* ^ H* + g J2r ^(-R); where $ is a scaling operator 
with scaling dimension x$ = 2. We may develop the current-current correla- 



tion function (|T8| in a power series in g, the coefficient of each term being a 
sum over the Rj of correlation functions ( J^(ri) Ji.(r2)$(-Ri)), ( J^(ri) J,^(r2)$(-Ri)$(-R2)), 
and so on, each evaluated with respect to the fixed point hamiltonian. The 
form of the ri and r2-dependence of each of these correlation functions is com- 
pletely fixed by conformal invariance in two dimensions, so that they may be 
computed in a simple model. Choosing a gaussian theory with hamiltonian 
H* = I a conserved current J^j ~ 9^0, and the marginal opera- 

tor $ ~ (90)^, all the correlators may be evaluated using Wick's theorem. 
For the 0{g) correction, it turns out that the only non-zero components (in 
complex coordinates) have n = z, u = z, and vice versa. The form of the 
correlation function is 

{Jz{zi)Mz2M0)) oc l/{zjzl) (52) 

where we have set i?i = for convenience. 

Now the 0{g) correction to the total area within all loops (|15D is 

g J {Jy{zi)Jy{z2)^{Ri))\xi - X2\6iyi - y2)dxidx2dyidy2(fRi (53) 

where Jy oc Jz — Jz- This is to be evaluated in a large but finite region of 
linear size 0{L). As before, we shall use the infinite volume continuum limit 
form (^) of the correlation function, justifying this a posteriori. The integral 
in (|53|) is then proportional to the area A of the system, and we remove this 
factor by setting Ri = 0. The remaining integral is then proportional to 

/ 7 , — dxidx2dyi (54) 

J-oo {xi +tyiy{x2 -tyiY 

The contour integration over yi vanishes unless Xi and X2 have the same sign: 
the result is then proportional to 



/ / 1—^ — ^^^^dxidx2= [^^r^ln{L/a) (55) 
Jo Jo ixi — X2)'* J Xi 



with an equal contribution from xi,X2 < 0. We have cut off the logarithmi- 
cally divergent integral in the last step, arguing that because the divergence 
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is only logarithmic, it was permissible to use the infinite-volume forms for 
the correlation function in the integrand. 

The important point about this result is that it is 0{glnL), not 0{g{lnLy ) 
as might have been expected (recall that the leading term is 0{g^lnL)). A 
similar, but more tedious, calculation shows that the next term is 0{g^\nL), 
and we conjecture that the nth order term is 0{g"' In L). This is consistent 
with the fact that, in the gaussian model, g is exactly marginal so that k, 
the coefficient of the O(lnL) term, depends continuously on g. 

However, in the 4-state Potts model the perturbation is not exactly 
marginal, and instead fiows logarithmically slowly to zero under the RG. 
This may be taken into account by replacing the bare expansion param- 
eter g by the running coupling 

g{L) = ^^-^ - {b\nL)-' + 0{g-\\nL)-') (56) 

where 6 is a known constant whose value is not important. 

Inserting this result into the formula (|T3p for the total area (Atot) ~ 
AEA<o{L^)NiA) gives 



E iV(^)~2ClnL(l + ^+|J|, + 0((l„L)- 



3\ 



(57) 



A<0{L'2) 



where the aj are non-universal constants. Differentiating this with respect 
to ~ A then gives the main result quoted in (|49|) 

^(^)~^(l-(l^ + 0((lnA)-)) (58) 

The interesting feature of this result is the absence of the 0{{\nA)~^) term, 
proportional to Oi. 
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